/*
 * GeoGeoReactor.cpp
 *
 *  Created on: 18 Oct 2011
 *      Author: allan
 */

#include "GeoReactor.h"

// Eigen includes
#include <Eigen/Core>
using namespace Eigen;

// GeoReact includes
#include "Geochemistry/GeoState.h"
#include "Geochemistry/GeoSystem.h"


GeoReactor::GeoReactor(const GeoSystem& system) :
system(system), dtMin(1.0E-6), dtMax(unsigned(-1)), dt(1.0E-6)
{
	kSolver.Initialiaze(system.GetMultiphase(), 
	                    system.GetEquilibriumReactions(), 
	                    system.GetKineticReactions(), 
	                    system.GetPrimarySpecies(), 
	                    system.GetEquilibriumSpecies(), 
	                    system.GetKineticSpecies());
	
	kSolver.SetControlSolution(0, 1.0E-6);
	
	kSolverEuler.Initialiaze(system.GetMultiphase(), 
	                         system.GetEquilibriumReactions(), 
	                         system.GetKineticReactions(), 
	                         system.GetPrimarySpecies(), 
	                         system.GetEquilibriumSpecies(), 
	                         system.GetKineticSpecies());
}

void GeoReactor::SetMinimumTimeStep(double dtMin)
{
	this->dtMin = dtMin;
	
	if(dt < dtMin) dt = dtMin;
}

void GeoReactor::SetMaximumTimeStep(double dtMax)
{
	this->dtMax = dtMax;
	
	if(dt > dtMax) dt = dtMax;
}

void GeoReactor::SetAlgorithm(OdeAlgorithm algorithm)
{
	kSolver.SetAlgorithm(algorithm);
}

void GeoReactor::SetAdaptiveTimeStepControl(double epsAbs, double epsRel)
{
	kSolver.SetControlSolution(epsAbs, epsRel);
}

const double GeoReactor::GetTimeStep() const
{
	return dt;
}

void GeoReactor::React(double& t, double tEnd, GeoState& state)
{
	  double T = state.GetTemperature();
	  double P = state.GetPressure();
	VectorXd n = state.GetComposition();
	
	kSolver.Iterate(t, tEnd, dt, T, P, n);
	
	if(dt < dtMin) dt = dtMin;
	if(dt > dtMax) dt = dtMax;
	
	state.SetComposition(n);
}

void GeoReactor::ReactEuler(double& t, double tEnd, GeoState& state)
{
	  double T = state.GetTemperature();
	  double P = state.GetPressure();
	VectorXd n = state.GetComposition();
	
	kSolverEuler.Iterate(t, tEnd, dt, T, P, n);
	
	if(dt < dtMin) dt = dtMin;
	if(dt > dtMax) dt = dtMax;
	
	state.SetComposition(n);
}
